############################################
############################################
# previous: 1-VAP_data #####################
############################################

###### Summary ############################################################################################
#This script extracts the data for each individual species in the master list from the full data file based on
#species names and known synonyms. It cross references the data to a data vetting spreadsheet that lists which 
#records or locations should be excluded for each species and which records need to be corrected from one species 
#to another species. It then removes the 25% least accurate unique occurrence records per 1km grid cell unless 
#the species is data deficient, in which case all data is used. The result is a number of files for each species, 
#including one with the original data, one with the cleaned data, one with the remaining data after removal of 
#inaccurate records, and one final file with only species, latitude and longitude for use with Maxent. Lastly,
#the script makes maps of remaining records for each species, as well as for each group of species that will 
#need to be modelled together (see methods for 3-VAP-MaxentPrep), and each higher level species complex or group.
###########################################################################################################

############################################
# next: 3-VAP_MaxentPrep ###################
############################################

##############################################################################################
#convert GBIF files from tab delimited into comma delimed files before using this script

#load required packages:
library(R.utils)
library(raster)
library(maptools)
library(rgdal)
library(stringr)
library(sf)

#set project name:
projectName<-"WHOSnakes"

#define directories:
MyDir<-paste("/home/jc217070/Maxent",sep="")
#MyDir<-paste("C:/Users/jc217070/Dropbox/temp/WHOSnakeData",sep="")
PlotDir<-paste(MyDir,"/Maxent_Data/Plots/Plots_WHOSnakes_2020_09_01",sep="")

#read species reference list:





SppSum<-read.table(file=file.path(paste(MyDir, "/Maxent_LUTables/",projectName,".csv",sep="")), header=TRUE, sep=",") #load data summary for all species
#SppSum<-read.table(file=file.path(paste(MyDir, "/Maxent_LUTables/",projectName,"_Fin1.csv",sep="")), header=TRUE, sep=",") #load data summary for all species







spp<-SppSum$gen_sp_subtax

Vet<-read.table(file=file.path(paste(MyDir, "/Maxent_LUTables/",projectName,"_Exclusions_Final.csv",sep="")), header=TRUE, sep=",") 
Vet$sp<-gsub(" ","_",Vet$sp)
Vet$shouldbe<-gsub(" ","_",Vet$shouldbe)

Vet2<-SppSum
spExcept<-Vet2$gen_sp_subtax[Vet2$spExcept==1]

#get spatial data layers for plots:
mybg<-raster(paste(PlotDir,"/Cop_FAPAR_mean_2010-2019_WW1km.asc",sep=""))
worldmap <- st_read(paste(PlotDir,"/worldmap/WHO_countries.shp",sep=""))


WHOpolies <- st_read(paste(MyDir,"/WHO snake distributions shapefiles Online Database 2021/MasterDistributionV2.shp",sep=""))
GARDpolies <- st_read(paste(MyDir,"/GARD ranges/modeled_reptiles.shp",sep=""))
WHOpoliesNEW <-st_read(paste(MyDir,"/WHOSnakesShapes/SSN_CTY_CAT_BRG_FEA_IN_Mar2023.shp",sep=""))
  
#mycols1<-c("cornsilk","khaki1","darkolivegreen1","olivedrab1","greenyellow",
#           "lawngreen","green2","green3","green4","darkgreen")
mycols1<-c("grey90","grey85","grey80","grey75","grey70",
           "grey65","grey60","grey55","grey50","grey45",
           "grey40","grey35","grey30")
mycols2<-c("royalblue","cyan","navy","skyblue","gold",
           "red","darkred","darkorange","goldenrod1",
           "deeppink","cornflowerblue","purple","slateblue1","lightsalmon","burlywood4","coral3",
           "deeppink4","darkorange4","cyan4","purple4")
WHOcol <- rgb(150, 180, 250, max = 255, alpha = 120)
WHOcol2 <- rgb(0, 150, 255, max = 255, alpha = 0)
WHOcol3 <- rgb(169, 160, 160, max = 255, alpha = 100)
WHOcol4 <- rgb(211,211,211, max = 255, alpha = 100)

#start of 'justspp' script
################################################################
################################################################
################################################################
# make first lot of files for each actual species separately ###
################################################################

alldata<-read.table(file=file.path(paste(MyDir, "/Maxent_Data/CLEAN_DATA/",projectName,"_allRecs.csv",sep="")), header=TRUE, sep=",") #load data summary for all species
ExtraData<-read.table(file=file.path(paste(MyDir, "/Maxent_Data/CLEAN_DATA/",projectName,"_georef_manually.csv",sep="")), header=TRUE, sep=",") #load data summary for all species

#make column for reference number for the category of record quality & origin: 1=not modified from original datasets, 
#2=verified by recent taxonomic study, 3=manually modified taxon ID or location
alldata$VetCom<-1 
alldata$VetCom<-ifelse(alldata$sourcefile=="LIT" & alldata$issue=="verified",2,1)

alldata$OUN<-ifelse(is.na(alldata$OUN) | alldata$OUN=="" , alldata$UN, alldata$OUN)

#make unique location columns to match record exclusions and additions to: rounded to 10 km,
#so anything within 10km of an excluded record or a similar addition record get found
alldata$unloc<-paste(round(alldata$lat,1),";",round(alldata$long,1),sep="")
#alldata$fulldeg<-paste(gsub("-","",as.character(paste(round(alldata$lat)))),";",gsub("-","",as.character(paste(round(alldata$long),sep=""))),sep="")
Vet$unloc<-paste(round(Vet$lat,1),";",round(Vet$long,1),sep="")

spp<-SppSum$gen_sp_subtax




###
sppB<-SppSum$gen_sp_subtax[SppSum$Redo==1]
#sppB<-c("Walterinnesia_aegyptia","Walterinnesia_morgani")
#sppB<-spp
###





allsynonyms<-vector()

for (sp in spp){
  n<-which(spp==sp)
  spsynonyms<-str_split(SppSum$Data_Synonyms[n], ";;")          #split START_DATE in 3 separate elements
  spsynonyms<-c(spsynonyms[[1]], paste(gsub(" 01","",SppSum$WHOSpecies[n])), paste(SppSum$Accepted_Species[n]),gsub("_"," ",sp))
  spsynonyms<-spsynonyms[spsynonyms!="NA"]
  allsynonyms<-c(allsynonyms, spsynonyms)
}







for (sp in spp){
if (sp %in% sppB){
  
#if (sp %in% c("Trimeresurus_albolabris","Trimeresurus_guoi","Trimeresurus_aff._albolabris","Trimeresurus_caudornatus","Trimeresurus_insularis")){
  
  
  
  
  
  
  
  n<-which(spp==sp)
  
  if(n==1 & length(spp)==length(sppB)){ 
    #make empty columns in summary file to paste species stats into:
    #(only reset all these if starting a full new run and 
    #if not just rerunning & replacing species number 1)
    SppSum$year<-as.character(paste("NA"))
    SppSum$precision<-as.character(paste("NA"))
    SppSum$rawrecs<-as.character(paste("NA"))
    SppSum$finrecs<-as.character(paste("NA"))
    
    SppSum$SPswdrecs<-as.character(paste("NA"))
    SppSum$maprecs<-as.character(paste("NA"))
    SppSum$MAPswdrecs<-as.character(paste("NA"))
    
  }
  
  myMapSpecies<-SppSum$ModUnit[n]
  spelems<-str_split_fixed(SppSum$WHOSpecies[n], " ", 3)
  WHOSpecies<-paste(spelems[,1],spelems[,2])
  ListSpecies<-SppSum$Listing_Name[n]
    
    
  Sys.sleep(0.1)
  print(paste(Sys.time()," extracting data for ", sp, ": ", n, " out of ", length(spp), sep=""))
  
  Sys.sleep(0.1)
  print(paste(Sys.time()," correct synonyms", sep=""))
  
  #create colum with corrected taxonomic name in full dataset before extracting data for this species
  synonyms<-str_split(SppSum$Data_Synonyms[n], ";;")          #split START_DATE in 3 separate elements
  synonyms<-c(synonyms[[1]], paste(gsub(" 01","",SppSum$WHOSpecies[n])), paste(SppSum$Accepted_Species[n]),gsub("_"," ",sp))
  synonyms<-synonyms[synonyms!="NA"]
  notsynonyms<-allsynonyms[!(allsynonyms %in% synonyms)]
  
  alldata$myTaxon<-ifelse((alldata$species %in% synonyms | alldata$scientificName %in% synonyms 
                           | paste(alldata$genus," ",alldata$species,sep="") %in% synonyms
                           | paste(alldata$genus," ",alldata$species," ",alldata$subtaxon, sep="") %in% synonyms)
                          & !(alldata$scientificName %in% notsynonyms),
                          paste(sp),paste(alldata$myTaxon))
  
  alldata$ModUnit<-ifelse(alldata$myTaxon==sp,paste(myMapSpecies),paste(alldata$ModUnit))
  
  Sys.sleep(0.1)
  print(paste(Sys.time()," apply corrections and exclusions", sep=""))
  
  Extrarecs<- unique(ExtraData[ExtraData$myTaxon==paste(sp),])
  sprecs_georef<- unique(ExtraData[ExtraData$myTaxon==paste(sp),])
  write.table(x=sprecs_georef, file=paste(MyDir, "/Maxent_SWDs/spp_occs/",sp,"_georef.csv",sep=""),qmethod = "double",
              na = "NA", row.names = FALSE, col.names = TRUE, sep=",",quote=TRUE)
  
  sprecso<- unique(alldata[alldata$myTaxon==paste(sp),])
  sprecso<-subset(sprecso, !is.na(sprecso$lat))
  sprecso<-subset(sprecso, !is.na(sprecso$long))
  
  if(sp== "Vipera_ammodytes"){
    sprecso<-subset(sprecso, sprecso$lat < 49)
  }
  
  
  
  
  
  
  
  if(length(sprecso[,1])>=1){
    
    write.table(x=sprecso, file=paste(MyDir, "/Maxent_SWDs/spp_occs/",sp,"_allRecsraw.csv",sep=""),qmethod = "double",
                na = "NA", row.names = FALSE, col.names = TRUE, sep=",",quote=TRUE)
  }
    VetExc<-Vet[Vet$sp==paste(sp)|Vet$sp==paste(gsub(" ","_",gsub(" 01","",SppSum$WHOSpecies[n]))),]       #records that need to be removed
    VetExc<-subset(VetExc,!(is.na(VetExc$sp)))
    
    VetAdd<-Vet[Vet$shouldbe==paste(sp)|Vet$shouldbe==paste(gsub(" ","_",gsub(" 01","",SppSum$WHOSpecies[n]))),] #records that need to be added
    VetAdd<-subset(VetAdd,!(is.na(VetAdd$shouldbe)))
    
    
    
    
    #exclude faulty records:
    sprecs<-sprecso[!(sprecso$OUN %in% VetExc$recID) & !(sprecso$OUN %in% VetExc$finalUN),] #first exclude records whos OUN is explicitly listed in my exclusions
    sprecsexcl1<-sprecso[(sprecso$OUN %in% VetExc$recID) | (sprecso$OUN %in% VetExc$finalUN),]
    
    
    
    
    
    sprecs<-sprecs[!(sprecs$unloc %in% VetExc$unloc),] #then exclude any records with the same locations as the excluded ones
    sprecsexcl2<-sprecs[(sprecs$unloc %in% VetExc$unloc),]
    
    sprecsexclall<-rbind(sprecsexcl1,sprecsexcl2)
    write.table(x=sprecsexclall, file=paste(MyDir, "/Maxent_SWDs/spp_occs/",sp,"_exclude.csv",sep=""),qmethod = "double",
                na = "NA", row.names = FALSE, col.names = TRUE, sep=",",quote=TRUE)

    #also get records from complete data that were the wrong species before and need to be switched over to this one
    #also get records from complete data that were the wrong species before and need to be switched over to this one
    if(length(VetAdd$sp)>0){
      
      #gett extra records from alldata
      sprecsExtra<-sprecs[1,] #placeholder line to bind rest to
      
      for (addsp in unique(VetAdd$sp)){ #for each species for which records need to be changed to my species
        
        VetAddsub<-VetAdd[VetAdd$sp==addsp,] #get the subset of additions that used to be this incorrect species
        #VetAddsub$tempUN<-seq(1:length(VetAddsub[,1]))
        
        if(addsp %in% SppSum$gen_sp_subtax){
          i<-which(SppSum$gen_sp_subtax==addsp)
          sprecsExtra2<-alldata[alldata$myTaxon == addsp & (alldata$OUN %in% unique(VetAddsub$recID)| alldata$OUN %in% unique(VetAddsub$finalUN)|alldata$unloc %in% VetAddsub$unloc),] #get records for this incorrect species from alldata
          sprecsExtra<-rbind(sprecsExtra,sprecsExtra2) #add these new records to the list of records that need to be added to my species
          
        }else{
          if(addsp %in% gsub(" ","_",gsub(" 01","",SppSum$WHOSpecies))){
            i<-which(gsub(" ","_",SppSum$WHOSpecies)==addsp)
            sprecsExtra2<-alldata[(alldata$myTaxon == addsp | alldata$myTaxon == SppSum$gen_sp_subtax[i]) & (alldata$OUN %in% unique(VetAddsub$recID) | alldata$OUN %in% unique(VetAddsub$finalUN)|alldata$unloc %in% VetAddsub$unloc),] #get records for this incorect species from alldata
            sprecsExtra<-rbind(sprecsExtra,sprecsExtra2) #add these new records to the list of records that need to be added to my species
          }else{
            sprecsExtra<-sprecsExtra #add these new records to the list of records that need to be added to my species
            print(paste("cannot find species fro which to transfer records over (",addsp,")!!!"))
          }}}
      
      #get rid of any double ups of records I want to add
      sprecsExtra<-sprecsExtra[2:length(sprecsExtra[,1]),] #get rid of first placeholder line
      sprecsExtra<-unique(sprecsExtra) #get rid of doubled up lines
      sprecsExtra<-subset(sprecsExtra, !is.na(sprecsExtra$lat)) #get rid of lines with no location data
      
      if(length(sprecsExtra[,1])>0){ #only do this if there's any records to be added
        
        #substitute faulty lat longs with correct lat longs from vetting file
        for (rec in c(1:length(sprecsExtra[,1]))){
          
          l<-which(VetAdd$unloc==sprecsExtra$unloc[rec] & VetAdd$sp==sprecsExtra$myTaxon & !(is.na(VetAdd$shouldlat)))
          #k<-which(VetAdd$recID==paste(sprecsExtra$OUN[rec]))
          
          if(length(l)>0){
            if(!(is.na(as.numeric(VetAdd$shouldlat[l][1])))){
              sprecsExtra$lat[rec]<-unique(VetAdd$shouldlat[l])[1]
              sprecsExtra$long[rec]<-unique(VetAdd$shouldlong[l])[1]
            }
          }
          
        }
        
        sprecsExtra$VetCom<-3
        sprecsExtra$myTaxon<-paste(sp)
        sprecsExtra$ModUnit<-paste(myMapSpecies)
        #update unloc
        sprecsExtra$unloc<-paste(round(sprecsExtra$lat,1),";",round(sprecsExtra$long,1),sep="")
        
        #bind extra data to file for this species
        sprecs<-rbind(sprecs,sprecsExtra)
        
      }
    }
    
    
    #exclude faulty records (needs to be repeated after additions in case records were 
    #added and later excluded; however only do for updated unloc because some excluded 
    #records were supposed to be added back in with new lat/long: if these were excluded 
    #again later, they would be excluded under the updated unloc and would still be found 
    sprecs<-sprecs[!(sprecs$OUN %in% VetExc$recID) & !(sprecs$OUN %in% VetExc$finalUN),] #first exclude records whos OUN is explicitly listed in my exclusions
    sprecs<-sprecs[!(sprecs$unloc %in% VetExc$unloc),] #then exclude any records with the same locations as the excluded ones
   
    Sys.sleep(0.1)
    print(paste(Sys.time()," reduce to final records", sep=""))
    
    #make all NA uncertainties 100km (not great but also not horrible) so that I can exclude worst records based on stats
    sprecs<- sprecs[order(-(sprecs$lattrunc1km), sprecs$longtrunc1km, (sprecs$uncertainty), sprecs$year), ] #sort by lat, long and moxt precise and most recent rows first
    sprecs$uncertainty<-ifelse(is.na(sprecs$uncertainty),100000,sprecs$uncertainty)
    maxunc<-max(na.omit(sprecs$uncertainty))
    minyear<-min(na.omit(sprecs$year))
    sprecs<-subset(sprecs, !is.na(sprecs$myTaxon))
    sprecs<-subset(sprecs, !is.na(sprecs$lat))
    sprecs<- sprecs[order(-(sprecs$lattrunc1km), sprecs$longtrunc1km, (sprecs$uncertainty), sprecs$year), ] #sort by lat, long and moxt precise and most recent rows first
    
    if (length(sprecs[,1])>0){
    
      unrecs<-sprecs[!duplicated(sprecs[,c("longtrunc1km", "lattrunc1km")]),]
      myQ75<-quantile(unrecs$uncertainty,probs = 0.75,na.rm=TRUE)
      n1<-length(unrecs[,1])
      n2<-length(unrecs[unrecs$uncertainty<=myQ75,1])
      
    #write out copy of species records before modifications:
    write.table(x=sprecs, file=paste(MyDir, "/Maxent_SWDs/spp_occs/",sp,"_allRecs.csv",sep=""),qmethod = "double",
                na = "NA", row.names = FALSE, col.names = TRUE, sep=",",quote=TRUE)

    
    if(sp %in% spExcept |is.na(myQ75)|n2<=20){ 
      #if the species is listed as an exclusion because of eg spatial bias in preision, use all records
      redrecs<-sprecs
      SppSum$year[n]<-paste(minyear)
      SppSum$precision[n]<-paste(maxunc)
      SppSum$rawrecs[n]<-paste(as.character(n1))
      SppSum$finrecs[n]<-paste(as.character(n1))
    } else {
      #exclude records by 75% quantile of precision (of unique records)
      redrecs<-sprecs[sprecs$uncertainty<=myQ75,]
      SppSum$year[n]<-paste(minyear)
      SppSum$precision[n]<-paste(myQ75)
      SppSum$rawrecs[n]<-paste(as.character(n1))
      SppSum$finrecs[n]<-paste(as.character(n2))
    }
    
    
    #only keep most precise and most recent (see sorting above) unique records:
    redrecs<-redrecs[!duplicated(redrecs[,c("longtrunc1km", "lattrunc1km")]),]
    
    #write out copy of reduced species records:
    redrecs<-subset(redrecs, !is.na(redrecs$myTaxon))
    redrecs<-subset(redrecs, !is.na(redrecs$lat))
    write.table(x=redrecs, file=paste(MyDir, "/Maxent_SWDs/spp_occs/",sp,"_redRecs.csv",sep=""),qmethod = "double",
                na = "NA", row.names = FALSE, col.names = TRUE, sep=",",quote=TRUE)
    
    #make final csv file with only taxon, lat & long:
    finrecs<-redrecs[,c("myTaxon","long", "lat")]
    
    #write out copy of reduced species records:
    finrecs<-subset(finrecs, !is.na(finrecs$myTaxon))
    finrecs<-subset(finrecs, !is.na(finrecs$lat))
    write.table(x=finrecs, file=paste(MyDir, "/Maxent_SWDs/spp_occs/",sp,".csv",sep=""),qmethod = "double",
                na = "NA", row.names = FALSE, col.names = TRUE, sep=",",quote=TRUE)
    
    #write out copy of new species summary file and read it back in:
    write.table(x=SppSum, file=paste(MyDir, "/Maxent_LUTables/",projectName,"_Fin1b.csv",sep=""),
                row.names = FALSE, col.names = TRUE, sep=",")
    SppSum<-read.table(file=file.path(paste(MyDir, "/Maxent_LUTables/",projectName,"_Fin1b.csv",sep="")), header=TRUE, sep=",") #load data summary for all species
    
    Sys.sleep(0.1)
    print(paste(Sys.time()," prep mapping environment", sep=""))
    
    #make extent object describing a box around occurrence records (+-5degrees):
    species_shapeA <- subset(WHOpolies, WHOpolies$NAME_SSN == gsub("_"," ",paste(sp)))
    species_shapeB <- subset(WHOpolies, WHOpolies$NAME_SSN == gsub("_"," ",paste(WHOSpecies)))
    species_shapeC <- subset(GARDpolies,GARDpolies$Binomial == gsub("_"," ",paste(sp)))
    species_shapeD <- subset(GARDpolies, GARDpolies$Binomial == gsub("_"," ",paste(WHOSpecies)))
    species_shapeE <- subset(WHOpoliesNEW, WHOpoliesNEW$Scientific == gsub("_"," ",paste(sp)))
    species_shapeF <- subset(WHOpoliesNEW, WHOpoliesNEW$Scientific == gsub("_"," ",paste(WHOSpecies)))
    
    #make extent object describing a box around occurrence records (margin= 4 times x and y extent of occurrence records):
    
    
    
    if (length(finrecs[,1])>=1){
      
    xrange<- (max(na.omit(finrecs$long)))-(min(na.omit(finrecs$long)))
    yrange<- (max(na.omit(finrecs$lat)))-(min(na.omit(finrecs$lat)))
    xcenter<- (min(na.omit(finrecs$long)))+(xrange/2)
    ycenter<- (min(na.omit(finrecs$lat)))+(yrange/2)
    rangemax<-max(xrange,yrange)
    myext<-extent(ifelse((xcenter-(rangemax/2)- 5) >= -180, (xcenter-(rangemax/2)- 5), -180),
                  ifelse((xcenter+(rangemax/2)+ 5) <=  180, (xcenter+(rangemax/2)+ 5),  180),
                  ifelse((ycenter-(rangemax/2)- 5) >=  -90, (ycenter-(rangemax/2)- 5),  -90),
                  ifelse((ycenter+(rangemax/2)+ 5) <=   90, (ycenter+(rangemax/2)+ 5),   90))
    } else{
      
      
      if(length(species_shapeE$OBJECTID)>=1){
        myext<-extent(species_shapeE)
        myext<-extent(myext[1]-5,myext[2]+5,myext[3]-5,myext[4]+5)
      }else{
        if(length(species_shapeF$OBJECTID)>=1){
          myext<-extent(species_shapeF)
          myext<-extent(myext[1]-5,myext[2]+5,myext[3]-5,myext[4]+5)
        } else{
          if(length(species_shapeA$OBJECTID)>=1){
            myext<-extent(species_shapeA)
            myext<-extent(myext[1]-5,myext[2]+5,myext[3]-5,myext[4]+5)
          }else{
            if(length(species_shapeB$OBJECTID)>=1){
              myext<-extent(species_shapeB)
              myext<-extent(myext[1]-5,myext[2]+5,myext[3]-5,myext[4]+5)
            } else{
              myext<-extent(min(na.omit(finrecs$long))-5,  
                            max(na.omit(finrecs$long))+5, 
                            min(na.omit(finrecs$lat))-5, 
                            max(na.omit(finrecs$lat))+5)
            }
          }
        }
      }
 
    }
    
    mybgxs<-crop(mybg, myext)

    i<-which(Vet2$gen_sp_subtax==sp)
    yesCountry<-str_split(Vet2[i,"country"], "/")                   #split countries of occurrence into separate elements
    mycountries <- subset(worldmap, worldmap$GU_A3 %in% yesCountry[[1]])     #get subset of countries that species occurs in

    Sys.sleep(0.1)
    print(paste(Sys.time()," make uncertainty buffers", sep=""))
    
    #make uncertainty buffers
    
    
    sprecs$uncertainty<-round((sprecs$uncertainty/10000),0)*10000
    sprecs$uncertainty<-ifelse(sprecs$uncertainty<=5000,5000,sprecs$uncertainty)
    sprecs$uncertainty<-ifelse(sprecs$uncertainty>=1000000,1000000,sprecs$uncertainty)
    
    w <- unique(sprecs$uncertainty)
    w<- sort(w) #sort small to large value
    b <- list()
    for (i in 1:length(w)) {
      x <- SpatialPoints(sprecs[sprecs$uncertainty==w[i],c(14,13)])
      b[[i]] <- buffer(x, width=w[i], dissolve=TRUE)
    } 
    if(length(b)>1){
      bb <- do.call(bind, b)
    } else{
      bb<-b[[1]]
    }
    
    bb2<-st_as_sf(bb)
    st_write(bb2, paste(PlotDir,"/Buffers/",sp,".shp",sep=""),delete_layer=TRUE)
    
    Sys.sleep(0.1)
    print(paste(Sys.time()," make maps", sep=""))
    
    #make PNG plots:
    png(filename=paste(PlotDir,"/",sp,"_recs.png",sep=""), units="in", width=8, height=12, res=500) #png smaller file size, A3 for better resolution
    par(mfrow = c(2,1), mar = c(2.5,2.5,0,0), oma=c(1,1,4,1), mgp=c(1.2,0.5,0))
    #mar: order: bottom, left, top, and right
    
    ################
    #1. Plot########
    ################
    
    #plot background
    plot(mybg, col=mycols1,box=FALSE, bty="n", cex.axis=0.8, legend=FALSE, axes=FALSE)
    plot(worldmap[5], border="grey20", lwd= 0.5, add=TRUE,col=NA,reset=FALSE)
    
    #plot old WHO distribution shape file
    if(!(length(species_shapeE$OBJECTID)>=1) & !(length(species_shapeF$OBJECTID)>=1)){ #if there's no updated polygons, plot old ones in pink
      
      if(length(species_shapeA$OBJECTID)>=1){
        plot(species_shapeA[2], add=TRUE, border="cornflowerblue",col=WHOcol,lwd= 1)
      }else{
        if(length(species_shapeB$OBJECTID)>=1){
          plot(species_shapeB[2], add=TRUE, border="cornflowerblue",col=WHOcol,lwd= 1)
        }
      }
      
      #plot GARD distribution shape file
      if(length(species_shapeC$Binomial)>=1){
        plot(species_shapeC[2], add=TRUE, border="red",col=WHOcol2,lwd= 1, lty=2)
      }else{
        if(length(species_shapeD$Binomial)>=1){
          plot(species_shapeD[2], add=TRUE, border="red",col=WHOcol2,lwd= 1, lty=2)
        }
      }
      
      
    } else { # if there ARE updated polygons available, plot old ones in grey, new ones in pink, and no GARD polies
      
      if(length(species_shapeA$OBJECTID)>=1){
      plot(species_shapeA[2], add=TRUE, border="darkgrey",col=WHOcol4,lwd= 1)
    }else{
      if(length(species_shapeB$OBJECTID)>=1){
        plot(species_shapeB[2], add=TRUE, border="darkgrey",col=WHOcol4,lwd= 1)
      }
    }
      #plot NEW WHO distribution shape file
      if(length(species_shapeE$OBJECTID)>=1){
        plot(species_shapeE[2], add=TRUE, border="cornflowerblue",col=WHOcol,lwd= 1)
      }else{
        if(length(species_shapeF$OBJECTID)>=1){
          plot(species_shapeF[2], add=TRUE, border="cornflowerblue",col=WHOcol,lwd= 1)
        }
      }
      
    }
    
    
    plot(mycountries[5], add=TRUE, col=NA, border="black",lwd= 1)

    #plot points
    points(sprecso$lat~sprecso$long, col="black", pch=1, cex=0.8)
    points(sprecs$lat~sprecs$long, col="darkred", pch=20, cex=0.9)
    points(sprecs$lat[sprecs$uncertainty<=100000]~sprecs$long[sprecs$uncertainty<=100000], col="orangered", pch=20, cex=0.9)
    points(sprecs$lat[sprecs$uncertainty<=20000]~sprecs$long[sprecs$uncertainty<=20000], col="royalblue", pch=20, cex=0.9)
    points(redrecs$lat~redrecs$long, col="black", pch=4, cex=1)
    #points(sprecs$lat[sprecs$VetCom==2]~sprecs$long[sprecs$VetCom==2], col="royalblue", pch=4, cex=1)
    
    text(-179,-95, labels=expression("__"),pos=4, col="cornflowerblue", cex=1.2)
    text(-164,-95.2, labels=expression("current WHO distribution estimate"),pos=4, col="black", cex=0.8)
    #text(-179,-105, labels=expression("__"),pos=4, col="darkgrey", cex=1.2)
    #text(-164,-105.2, labels=paste("old WHO dstribution estimate"),pos=4, col="black", cex=0.8)
    text(-179,-105, labels=expression("__"),pos=4, col="black", cex=1)
    text(-164,-105.2, labels=paste("countries where species MAY occur"),pos=4, col="black", cex=0.8)
    
    text(0,-95, labels=expression("."),pos=4, col="royalblue", cex=4)
    text(15,-95.2, labels=expression("precise records (<20km uncertainty)"),pos=4, col="black", cex=0.8)
    text(0,-105, labels=expression("."),pos=4, col="orangered", cex=4)
    text(15,-105.2, labels=paste("records with uncertainty 20-100km or NA"),pos=4, col="black", cex=0.8)
    text(0,-115, labels=expression("."),pos=4, col="darkred", cex=4)
    text(15,-115.2, labels=paste("high spatial uncertainty records"),pos=4, col="black", cex=0.8)
    text(4,-125.2, labels=expression("o"),pos=4, col="black", cex=0.8)
    text(15,-125.2, labels=paste("excluded records"),pos=4, col="black", cex=0.8)
    text(4,-135.2, labels=expression("x"),pos=4, col="black", cex=1)
    text(15,-135.2, labels=paste("Records used for SDM"),pos=4, col="black", cex=0.8)
    
    #################
    #2. Plot#########
    #################
    
    #plot background
    plot(mybgxs, col=mycols1,box=FALSE, bty="n", cex.axis=0.8, legend=FALSE, axes=FALSE)
    plot(worldmap[5], border="grey20", lwd= 0.5, add=TRUE,col=NA,reset=FALSE)
    
    if(!(length(species_shapeE$OBJECTID)>=1) & !(length(species_shapeF$OBJECTID)>=1)){ #if there's no updated polygons, plot old ones in pink
      
      if(length(species_shapeA$OBJECTID)>=1){
        plot(species_shapeA[2], add=TRUE, border="cornflowerblue",col=WHOcol,lwd= 1)
      }else{
        if(length(species_shapeB$OBJECTID)>=1){
          plot(species_shapeB[2], add=TRUE, border="cornflowerblue",col=WHOcol,lwd= 1)
        }
      }
      
      #plot GARD distribution shape file
      if(length(species_shapeC$Binomial)>=1){
        plot(species_shapeC[2], add=TRUE, border="red",col=WHOcol2,lwd= 1, lty=2)
      }else{
        if(length(species_shapeD$Binomial)>=1){
          plot(species_shapeD[2], add=TRUE, border="red",col=WHOcol2,lwd= 1, lty=2)
        }
      }
      
      
    } else { # if there ARE updated polygons available, plot old ones in grey, new ones in pink, and no GARD polies
      
      if(length(species_shapeA$OBJECTID)>=1){
        plot(species_shapeA[2], add=TRUE, border="darkgrey",col=WHOcol4,lwd= 1)
      }else{
        if(length(species_shapeB$OBJECTID)>=1){
          plot(species_shapeB[2], add=TRUE, border="darkgrey",col=WHOcol4,lwd= 1)
        }
      }
      #plot NEW WHO distribution shape file
      if(length(species_shapeE$OBJECTID)>=1){
        plot(species_shapeE[2], add=TRUE, border="cornflowerblue",col=WHOcol,lwd= 1)
      }else{
        if(length(species_shapeF$OBJECTID)>=1){
          plot(species_shapeF[2], add=TRUE, border="cornflowerblue",col=WHOcol,lwd= 1)
        }
      }
      
    }
    
    #plot(mycountries, add=TRUE, border="black",lwd= 1)
    plot(mycountries[5], add=TRUE, col=NA, border="black",lwd= 1)
    
    #plot(bb,border="grey60",col=WHOcol3, lwd=2,add=TRUE)
    
    #plot points
    points(sprecso$lat~sprecso$long, col="black", pch=1, cex=0.7)
    points(sprecs$lat~sprecs$long, col="darkred", pch=20, cex=0.9)
    points(sprecs$lat[sprecs$uncertainty<=100000]~sprecs$long[sprecs$uncertainty<=100000], col="orangered", pch=20, cex=0.9)
    #points(sprecs$lat[sprecs$uncertainty<=50000]~sprecs$long[sprecs$uncertainty<=50000], col="goldenrod1", pch=20, cex=0.9)
    points(sprecs$lat[sprecs$uncertainty<=20000]~sprecs$long[sprecs$uncertainty<=20000], col="royalblue", pch=20, cex=0.9)
    #points(sprecs$lat[sprecs$VetCom==2]~sprecs$long[sprecs$VetCom==2], col="royalblue", pch=3, cex=1)
    points(redrecs$lat~redrecs$long, col="black", pch=4, cex=1)
    
    #plot title in top margin of page:
    mtext(paste(gsub("_"," ",sp),sep=""), font= 4, side = 3, line = 1, outer = TRUE, cex= 2)
    mtext(paste("(WHO listed as '",ListSpecies,"')",sep=""), font= 3, side = 3, line = -0.4, outer = TRUE, cex= 1.5)

    #mtext(paste(SppSum$Comments_Taxonomy[n]), font= 1, side = 3, line = -1.9, outer = TRUE, cex= 1)
    #mtext(paste(SppSum$Comments_Range[n]), font= 1, side = 3, line = -2.8, outer = TRUE, cex= 1)
    #mtext(paste(SppSum$Comments_Records[n]), font= 1, side = 3, line = -3.7, outer = TRUE, cex= 1)
    #if(SppSum$Comments_Other[n]!="Other: no comment"){
    #mtext(paste(SppSum$Comments_Other[n]), font= 1, side = 3, line = -4.6, outer = TRUE, cex= 1)
    #}
    #The font "face" (1=plain, 2=bold, 3=italic, 4=bold-italic)
    
    dev.off()
    
  }
 }
#}
}













############################
# make second lot of plots for groups of species mapped together
############################

################
################
mapspp<-unique(SppSum$ModUnit)
#mapsppB<-unique(SppSum$ModUnit)
mapsppB<-unique(SppSum$ModUnit[SppSum$Redo==1])
#mapsppB<-c("Walterinnesia aegyptia","Walterinnesia morgani")
################
################

SppSum$maprecs<-paste("NA")

for (mapsp in mapspp){
  if(mapsp %in% mapsppB){
    
  ns<-which(SppSum$ModUnit==mapsp)
  m<-which(mapspp==mapsp)
  maprecs<-data.frame()
  mapname<-unique(SppSum$ModUnit)[m]
  
  Sys.sleep(0.1)
  print(paste(Sys.time()," plotting data for SDM modelling group ", mapsp, ": ", m, " out of ", length(mapspp), sep=""))
  
  #make PDF plots:
  png(filename=paste(PlotDir,"/SDMModGroup_",gsub(" ","_",mapsp),"_recs.png",sep=""), units="in", width=8.27, height=11.69, res=500) #png smaller file size, A3 for better resolution
  par(mfrow = c(2,1), mar = c(2.5,2.5,0,0), oma=c(1,1,4,3), mgp=c(1.2,0.5,0))
  
  plot(mybg, col=mycols1,box=FALSE, bty="n", cex.axis=0.8, legend=FALSE, axes=FALSE)
  plot(worldmap[5], border="grey20", lwd= 0.5, add=TRUE,col=NA,reset=FALSE)
  
  for (mapunit in SppSum[ns,"gen_sp_subtax"]){
    i<-which(SppSum[ns,"gen_sp_subtax"]==mapunit)
    k<-which(SppSum$gen_sp_subtax==mapunit)
    yesCountry<-str_split(SppSum[k,"country"], "/")                   #split countries of occurrence into separate elements
    mycountries <- subset(worldmap, worldmap$GU_A3 %in% yesCountry[[1]])     #get subset of countries that species occurs in
    plot(mycountries[5], add=TRUE, col=NA, border="black",lwd= 1)
  }
  
  #bind all records for species included in this map species together
  for (mapunit in SppSum[ns,"gen_sp_subtax"]){
    i<-which(SppSum[ns,"gen_sp_subtax"]==mapunit)
    
    if(file.exists(paste(MyDir, "/Maxent_SWDs/spp_occs/",mapunit,".csv",sep=""))){
      addrecs<-read.table(paste(MyDir, "/Maxent_SWDs/spp_occs/",mapunit,".csv",sep=""), sep=",", header=TRUE)
      maprecs<-rbind(maprecs,addrecs)
      
      redrecs<-read.table(paste(MyDir, "/Maxent_SWDs/spp_occs/",mapunit,"_redRecs.csv",sep=""), sep=",", header=TRUE)
      points(redrecs$lat~redrecs$long, col=mycols2[i], pch=20, cex=1)
      points(redrecs$lat[redrecs$VetCom==2]~redrecs$long[redrecs$VetCom==2], col=mycols2[i], pch=4, cex=1)
  
      #add a species label to the centre of the record cloud for each mapping unit
      myx<-mean(redrecs$long)
      myy<-mean(redrecs$lat)
      #graphics::text(x=myx, y=myy, labels = gsub("_"," ",mapunit), cex = 1, col = mycols2[i])
      
    }}
  
  if(length(maprecs)>0){
    l<-length(maprecs[,1])
    SppSum$maprecs[ns]<-paste(as.character(l))
    
    #write out copy of reduced species records:
    maprecs<-subset(maprecs, !is.na(maprecs$myTaxon))
    write.table(x=maprecs, file=paste(MyDir, "/Maxent_SWDs/spp_occs/",gsub(" ","_",mapsp),"_mod.csv",sep=""),qmethod = "double",
                na = "NA", row.names = FALSE, col.names = TRUE, sep=",",quote=TRUE)
    
    #make extent object describing a box around occurrence records (+-5degrees):
    myext<-extent(min(na.omit(maprecs$long))-5,  
                  max(na.omit(maprecs$long))+5, 
                  min(na.omit(maprecs$lat))-5, 
                  max(na.omit(maprecs$lat))+5)
    
  } else {
    myext<-extent(mybg)
  }
  
  #write out copy of new species summary file and read it back in:
  write.table(x=SppSum, file=paste(MyDir, "/Maxent_LUTables/",projectName,"_Fin1b.csv",sep=""),
              row.names = FALSE, col.names = TRUE, sep=",")
  SppSum<-read.table(file=file.path(paste(MyDir, "/Maxent_LUTables/",projectName,"_Fin1b.csv",sep="")), header=TRUE, sep=",") #load data summary for all species
  
  mybgxs<-crop(mybg, y=myext)
  plot(mybgxs, col=mycols1, box=FALSE, bty="n", cex.axis=0.8, legend=FALSE, axes=FALSE)
  plot(worldmap[5], border="grey20", lwd= 0.5, add=TRUE,col=NA,reset=FALSE)
  
  for (mapunit in SppSum[ns,"gen_sp_subtax"]){
    i<-which(SppSum[ns,"gen_sp_subtax"]==mapunit)
    k<-which(SppSum$gen_sp_subtax==mapunit)
    yesCountry<-str_split(SppSum[k,"country"], "/")                   #split countries of occurrence into separate elements
    mycountries <- subset(worldmap, worldmap$GU_A3 %in% yesCountry[[1]])     #get subset of countries that species occurs in
    plot(mycountries[5], add=TRUE, col=NA, border="black",lwd= 1)
  
    species_shapeA <- subset(WHOpolies, WHOpolies$NAME_SSN == gsub("_"," ",paste(mapunit)))
    #species_shapeB <- subset(WHOpolies, WHOpolies$NAME_SSN == gsub("_"," ",gsub(" 01","",paste(SppSum$WHOSpecies[k]))))
    species_shapeE <- subset(WHOpoliesNEW, WHOpoliesNEW$Scientific == gsub("_"," ",paste(mapunit)))
    #species_shapeF <- subset(WHOpoliesNEW, WHOpoliesNEW$Scientific == gsub("_"," ",gsub(" 01","",paste(SppSum$WHOSpecies[k]))))
    
      #if(!(length(species_shapeE$OBJECTID)>=1) & !(length(species_shapeF$OBJECTID)>=1)){ #if there's no updated polygons, plot old ones
      if(!(length(species_shapeE$OBJECTID)>=1)){ #if there's no updated polygons, plot old ones
        
        if(length(species_shapeA$OBJECTID)>=1){
          plot(species_shapeA[2], add=TRUE, border=mycols2[i],col=WHOcol,lwd= 1)
        }#else{
        #  if(length(species_shapeB$OBJECTID)>=1){
        #    plot(species_shapeB[2], add=TRUE, border=mycols2[i],col=WHOcol,lwd= 1)
        #  }
        #}
      } else { #otherwise use new versions
        
        if(length(species_shapeE$OBJECTID)>=1){
          plot(species_shapeE[2], add=TRUE, border=mycols2[i], col=WHOcol, lwd= 1)
        }#else{
        #if(length(species_shapeF$OBJECTID)>=1){
        #  plot(species_shapeF[2], add=TRUE, border=mycols2[i], col=WHOcol, lwd= 1)
        # }
        #}
      }}
    
  
  for (mapunit in SppSum[ns,"gen_sp_subtax"]){
    
    i<-which(SppSum[ns,"gen_sp_subtax"]==mapunit)
    if(file.exists(paste(MyDir, "/Maxent_SWDs/spp_occs/",mapunit,".csv",sep=""))){
      
      redrecs<-read.table(paste(MyDir, "/Maxent_SWDs/spp_occs/",mapunit,"_redRecs.csv",sep=""), sep=",", header=TRUE)
      points(redrecs$lat~redrecs$long, col=mycols2[i], pch=20, cex=1)
      points(redrecs$lat[redrecs$VetCom==2]~redrecs$long[redrecs$VetCom==2], col=mycols2[i], pch=4, cex=1)
      
      #add a species label to the centre of the record cloud for each mapping unit
      myx<-mean(redrecs$long)
      myy<-mean(redrecs$lat)
      graphics::text(x=myx, y=myy, labels = gsub("_"," ",mapunit), cex = 1, col = mycols2[i])
      
    }}
  
  #plot title in top margin of page:
  mtext(paste("SDM model group:"), font= 2, side = 3, line = 1, outer = TRUE, cex= 1.8)
  mtext(paste("'",mapname,"'",sep=""), font= 4, side = 3, line = -0.3, outer = TRUE, cex= 1.8)
  mtext(paste("includes species:"), font= 1, side = 3, line = -1.5, outer = TRUE, cex= 1)
  
  if(length(ns)<=4){
    mtext(paste(SppSum[ns,"Accepted_Species"], collapse = ', '), 
          font= 3, side = 3, line = -2.4, outer = TRUE, cex= 0.9)
  } else {
    if(length(ns)<=8){
    mtext(paste(SppSum[ns[1:4],"Accepted_Species"], collapse = ', '), 
          font= 3, side = 3, line = -2.4, outer = TRUE, cex= 0.9)
    mtext(paste(SppSum[ns[5:length(ns)],"Accepted_Species"], collapse = ', '), 
          font= 3, side = 3, line = -3.5, outer = TRUE, cex= 0.9)
    } else {
      if(length(ns)<=12){
        mtext(paste(SppSum[ns[1:4],"Accepted_Species"], collapse = ', '), 
              font= 3, side = 3, line = -2.4, outer = TRUE, cex= 0.9)
        mtext(paste(SppSum[ns[5:8],"Accepted_Species"], collapse = ', '), 
              font= 3, side = 3, line = -3.3, outer = TRUE, cex= 0.9)
        mtext(paste(SppSum[ns[9:length(ns)],"Accepted_Species"], collapse = ', '), 
              font= 3, side = 3, line = -4.2, outer = TRUE, cex= 0.9)
      } else {
        mtext(paste(SppSum[ns[1:4],"Accepted_Species"], collapse = ', '), 
              font= 3, side = 3, line = -2.4, outer = TRUE, cex= 0.9)
        mtext(paste(SppSum[ns[5:8],"Accepted_Species"], collapse = ', '), 
              font= 3, side = 3, line = -3.3, outer = TRUE, cex= 0.9)
        mtext(paste(SppSum[ns[9:12],"Accepted_Species"], collapse = ', '), 
              font= 3, side = 3, line = -4.2, outer = TRUE, cex= 0.9)
        mtext(paste(SppSum[ns[13:length(ns)],"Accepted_Species"], collapse = ', '), 
              font= 3, side = 3, line = -5.1, outer = TRUE, cex= 0.9)
    }}}
  #fonts: 1=normal, 2=bold, 3=cursive, 4= bold & cursive
  
  dev.off()
  
}
}









############################
# make third lot of overview plots for species complexes
############################

################
################
mapspp<-unique(SppSum$PlotGroup)
mapspp<-gsub(" ","_",gsub("; ","_",mapspp))
mapsppB<-unique(SppSum$PlotGroup[SppSum$Redo==1])
#mapsppB<-c("Walterinnesia spp.")
#mapsppB<-unique(SppSum$PlotGroup)
mapsppB<-gsub(" ","_",gsub("; ","_",mapsppB))

################
################

#SppSum$maprecs<-paste("NA")

for (mapsp in mapspp){
  if(mapsp %in% mapsppB){
    
  ns<-which(SppSum$PlotGroup==gsub("_"," ",mapsp))
  m<-which(mapspp==mapsp)
  maprecs<-data.frame()
  mapname<-unique(SppSum$PlotGroup)[m]
  
  Sys.sleep(0.1)
  print(paste(Sys.time()," plotting data for taxonomic species group ", gsub("_"," ",mapsp), ": ", m, " out of ", length(mapspp), sep=""))
  
  #make PDF plots:
  png(filename=paste(PlotDir,"/TaxGroup_",mapsp,"_recs.png",sep=""), units="in", width=8.27, height=11.69, res=500) #png smaller file size, A3 for better resolution
  par(mfrow = c(2,1), mar = c(2.5,2.5,0,0), oma=c(1,1,4,3), mgp=c(1.2,0.5,0))
  
  plot(mybg, col=mycols1,box=FALSE, bty="n", cex.axis=0.8, legend=FALSE, axes=FALSE)
  plot(worldmap[5], border="grey20", lwd= 0.5, add=TRUE,col=NA,reset=FALSE)
  
  for (mapunit in SppSum[ns,"gen_sp_subtax"]){
    i<-which(SppSum[ns,"gen_sp_subtax"]==mapunit)
    k<-which(SppSum$gen_sp_subtax==mapunit)
    yesCountry<-str_split(SppSum[k,"country"], "/")                   #split countries of occurrence into separate elements
    mycountries <- subset(worldmap, worldmap$GU_A3 %in% yesCountry[[1]])     #get subset of countries that species occurs in
    plot(mycountries[5], add=TRUE, col=NA, border="black",lwd= 1)
  }
  
  #bind all records for species included in this map species together
  for (mapunit in SppSum[ns,"gen_sp_subtax"]){
    i<-which(SppSum[ns,"gen_sp_subtax"]==mapunit)
    
    if(file.exists(paste(MyDir, "/Maxent_SWDs/spp_occs/",mapunit,".csv",sep=""))){
      addrecs<-read.table(paste(MyDir, "/Maxent_SWDs/spp_occs/",mapunit,".csv",sep=""), sep=",", header=TRUE)
      maprecs<-rbind(maprecs,addrecs)
      
      redrecs<-read.table(paste(MyDir, "/Maxent_SWDs/spp_occs/",mapunit,"_redRecs.csv",sep=""), sep=",", header=TRUE)
      points(redrecs$lat~redrecs$long, col=mycols2[i], pch=20, cex=1)
      points(redrecs$lat[redrecs$VetCom==2]~redrecs$long[redrecs$VetCom==2], col=mycols2[i], pch=4, cex=1)
      
      #add a species label to the centre of the record cloud for each mapping unit
      myx<-mean(redrecs$long)
      myy<-mean(redrecs$lat)
      #graphics::text(x=myx, y=myy, labels = gsub("_"," ",mapunit), cex = 1, col = mycols2[i])
      
    }}
  
  if(length(maprecs)>0){
    l<-length(maprecs[,1])
    #SppSum$maprecs[ns]<-paste(as.character(l))
    
    #write out copy of reduced species records:
    maprecs<-subset(maprecs, !is.na(maprecs$myTaxon))
    #write.table(x=maprecs, file=paste(MyDir, "/Maxent_SWDs/spp_occs/",mapsp,"_mod.csv",sep=""),qmethod = "double",
    #            na = "NA", row.names = FALSE, col.names = TRUE, sep=",",quote=TRUE)
    
    #make extent object describing a box around occurrence records (+-5degrees):
    myext<-extent(min(na.omit(maprecs$long))-5,  
                  max(na.omit(maprecs$long))+5, 
                  min(na.omit(maprecs$lat))-5, 
                  max(na.omit(maprecs$lat))+5)
    
  } else {
    myext<-extent(mybg)
  }
  
  mybgxs<-crop(mybg, y=myext)
  plot(mybgxs, col=mycols1, box=FALSE, bty="n", cex.axis=0.8, legend=FALSE, axes=FALSE)
  plot(worldmap[5], border="grey20", lwd= 0.5, add=TRUE,col=NA,reset=FALSE)
  
  for (mapunit in SppSum[ns,"gen_sp_subtax"]){
    i<-which(SppSum[ns,"gen_sp_subtax"]==mapunit)
    k<-which(SppSum$gen_sp_subtax==mapunit)
    yesCountry<-str_split(SppSum[k,"country"], "/")                   #split countries of occurrence into separate elements
    mycountries <- subset(worldmap, worldmap$GU_A3 %in% yesCountry[[1]])     #get subset of countries that species occurs in
    plot(mycountries[5], add=TRUE, col=NA, border="black",lwd= 1)
    
    species_shapeA <- subset(WHOpolies, WHOpolies$NAME_SSN == gsub("_"," ",paste(mapunit)))
    #species_shapeB <- subset(WHOpolies, WHOpolies$NAME_SSN == gsub("_"," ",gsub(" 01","",paste(SppSum$WHOSpecies[k]))))
    species_shapeE <- subset(WHOpoliesNEW, WHOpoliesNEW$Scientific == gsub("_"," ",paste(mapunit)))
    #species_shapeF <- subset(WHOpoliesNEW, WHOpoliesNEW$Scientific == gsub("_"," ",gsub(" 01","",paste(SppSum$WHOSpecies[k]))))
    
    #if(!(length(species_shapeE$OBJECTID)>=1) & !(length(species_shapeF$OBJECTID)>=1)){ #if there's no updated polygons, plot old ones
    if(!(length(species_shapeE$OBJECTID)>=1)){ #if there's no updated polygons, plot old ones
        
      if(length(species_shapeA$OBJECTID)>=1){
        plot(species_shapeA[2], add=TRUE, border=mycols2[i],col=WHOcol,lwd= 1)
      }#else{
       #  if(length(species_shapeB$OBJECTID)>=1){
       #    plot(species_shapeB[2], add=TRUE, border=mycols2[i],col=WHOcol,lwd= 1)
       #  }
       #}
    } else { #otherwise use new versions
      
      if(length(species_shapeE$OBJECTID)>=1){
        plot(species_shapeE[2], add=TRUE, border=mycols2[i], col=WHOcol, lwd= 1)
      }#else{
       #if(length(species_shapeF$OBJECTID)>=1){
       #  plot(species_shapeF[2], add=TRUE, border=mycols2[i], col=WHOcol, lwd= 1)
       # }
       #}
    }}
    
  for (mapunit in SppSum[ns,"gen_sp_subtax"]){
    
    i<-which(SppSum[ns,"gen_sp_subtax"]==mapunit)
    if(file.exists(paste(MyDir, "/Maxent_SWDs/spp_occs/",mapunit,".csv",sep=""))){
      
      redrecs<-read.table(paste(MyDir, "/Maxent_SWDs/spp_occs/",mapunit,"_redRecs.csv",sep=""), sep=",", header=TRUE)
      points(redrecs$lat~redrecs$long, col=mycols2[i], pch=20, cex=1)
      points(redrecs$lat[redrecs$VetCom==2]~redrecs$long[redrecs$VetCom==2], col=mycols2[i], pch=4, cex=1)
      
      #add a species label to the centre of the record cloud for each mapping unit
      myx<-mean(redrecs$long)
      myy<-mean(redrecs$lat)
      graphics::text(x=myx, y=myy, labels = gsub("_"," ",mapunit), cex = 1, col = mycols2[i])
      
    }}
  
  #plot title in top margin of page:
  mtext(paste("Taxonomic species group:"), font= 2, side = 3, line = 1, outer = TRUE, cex= 1.8)
  mtext(paste("'",mapname,"'",sep=""), font= 4, side = 3, line = -0.3, outer = TRUE, cex= 1.8)
  mtext(paste("includes species:"), font= 1, side = 3, line = -1.5, outer = TRUE, cex= 1)
  
  if(length(ns)<=4){
    mtext(paste(SppSum[ns,"Accepted_Species"], collapse = ', '), 
          font= 3, side = 3, line = -2.4, outer = TRUE, cex= 0.9)
  } else {
    if(length(ns)<=8){
      mtext(paste(SppSum[ns[1:4],"Accepted_Species"], collapse = ', '), 
            font= 3, side = 3, line = -2.4, outer = TRUE, cex= 0.9)
      mtext(paste(SppSum[ns[5:length(ns)],"Accepted_Species"], collapse = ', '), 
            font= 3, side = 3, line = -3.5, outer = TRUE, cex= 0.9)
    } else {
      if(length(ns)<=12){
        mtext(paste(SppSum[ns[1:4],"Accepted_Species"], collapse = ', '), 
              font= 3, side = 3, line = -2.4, outer = TRUE, cex= 0.9)
        mtext(paste(SppSum[ns[5:8],"Accepted_Species"], collapse = ', '), 
              font= 3, side = 3, line = -3.3, outer = TRUE, cex= 0.9)
        mtext(paste(SppSum[ns[9:length(ns)],"Accepted_Species"], collapse = ', '), 
              font= 3, side = 3, line = -4.2, outer = TRUE, cex= 0.9)
      } else {
        mtext(paste(SppSum[ns[1:4],"Accepted_Species"], collapse = ', '), 
              font= 3, side = 3, line = -2.4, outer = TRUE, cex= 0.9)
        mtext(paste(SppSum[ns[5:8],"Accepted_Species"], collapse = ', '), 
              font= 3, side = 3, line = -3.3, outer = TRUE, cex= 0.9)
        mtext(paste(SppSum[ns[9:12],"Accepted_Species"], collapse = ', '), 
              font= 3, side = 3, line = -4.2, outer = TRUE, cex= 0.9)
        mtext(paste(SppSum[ns[13:length(ns)],"Accepted_Species"], collapse = ', '), 
              font= 3, side = 3, line = -5.1, outer = TRUE, cex= 0.9)
      }}}
  #fonts: 1=normal, 2=bold, 3=cursive, 4= bold & cursive
  
  dev.off()
  
}
}







############################
# make third lot of overview plots for species merging
############################

################
################
#mapspp<-unique(SppSum$MergeSpecies[SppSum$PrimRegion=="Africa"])
mapspp<-unique(SppSum$MergeShort)
mapspp<-gsub(" ","_",gsub("; ","_",mapspp))
mapsppB<-unique(SppSum$MergeShort[SppSum$Redo==1])
#mapsppB<-c("Walterinnesia spp.")
#mapsppB<-unique(SppSum$MergeShort)
mapsppB<-gsub(" ","_",gsub("; ","_",mapsppB))
################
################

#SppSum$maprecs<-paste("NA")

for (mapsp in mapspp){
  if(mapsp %in% mapsppB){
  
  ns<-which(SppSum$MergeShort==gsub("_"," ",mapsp))
  m<-which(mapspp==mapsp)
  maprecs<-data.frame()
  
  mapname<-unique(SppSum$MergeSpecies)[m]
  
  Sys.sleep(0.1)
  print(paste(Sys.time()," plotting data for antivenom & combined impact group ", gsub("_"," ",mapsp), ": ", m, " out of ", length(mapspp), sep=""))
  
  #make PDF plots:
  png(filename=paste(PlotDir,"/AVImpactGroup_",mapsp,"_recs.png",sep=""), units="in", width=8.27, height=11.69, res=500) #png smaller file size, A3 for better resolution
  par(mfrow = c(2,1), mar = c(2.5,2.5,0,0), oma=c(1,1,4,3), mgp=c(1.2,0.5,0))
  
  plot(mybg, col=mycols1,box=FALSE, bty="n", cex.axis=0.8, legend=FALSE, axes=FALSE)
  plot(worldmap[5], border="grey20", lwd= 0.5, add=TRUE,col=NA,reset=FALSE)
  
  for (mapunit in SppSum[ns,"gen_sp_subtax"]){
    i<-which(SppSum[ns,"gen_sp_subtax"]==mapunit)
    k<-which(SppSum$gen_sp_subtax==mapunit)
    yesCountry<-str_split(SppSum[k,"country"], "/")                   #split countries of occurrence into separate elements
    mycountries <- subset(worldmap, worldmap$GU_A3 %in% yesCountry[[1]])     #get subset of countries that species occurs in
    plot(mycountries[5], add=TRUE, col=NA, border="black",lwd= 1)
  }
  
  #bind all records for species included in this map species together
  for (mapunit in SppSum[ns,"gen_sp_subtax"]){
    i<-which(SppSum[ns,"gen_sp_subtax"]==mapunit)
    
    if(file.exists(paste(MyDir, "/Maxent_SWDs/spp_occs/",mapunit,".csv",sep=""))){
      addrecs<-read.table(paste(MyDir, "/Maxent_SWDs/spp_occs/",mapunit,".csv",sep=""), sep=",", header=TRUE)
      maprecs<-rbind(maprecs,addrecs)
      
      redrecs<-read.table(paste(MyDir, "/Maxent_SWDs/spp_occs/",mapunit,"_redRecs.csv",sep=""), sep=",", header=TRUE)
      points(redrecs$lat~redrecs$long, col=mycols2[i], pch=20, cex=1)
      points(redrecs$lat[redrecs$VetCom==2]~redrecs$long[redrecs$VetCom==2], col=mycols2[i], pch=4, cex=1)
      
      #add a species label to the centre of the record cloud for each mapping unit
      myx<-mean(redrecs$long)
      myy<-mean(redrecs$lat)
      #graphics::text(x=myx, y=myy, labels = gsub("_"," ",mapunit), cex = 1, col = mycols2[i])
      
    }}
  
  if(length(maprecs)>0){
    l<-length(maprecs[,1])
    #SppSum$maprecs[ns]<-paste(as.character(l))
    
    #write out copy of reduced species records:
    maprecs<-subset(maprecs, !is.na(maprecs$myTaxon))
    #write.table(x=maprecs, file=paste(MyDir, "/Maxent_SWDs/spp_occs/",mapsp,"_mod.csv",sep=""),qmethod = "double",
    #            na = "NA", row.names = FALSE, col.names = TRUE, sep=",",quote=TRUE)
    
    #make extent object describing a box around occurrence records (+-5degrees):
    myext<-extent(min(na.omit(maprecs$long))-5,  
                  max(na.omit(maprecs$long))+5, 
                  min(na.omit(maprecs$lat))-5, 
                  max(na.omit(maprecs$lat))+5)
    
  } else {
    myext<-extent(mybg)
  }
  
  mybgxs<-crop(mybg, y=myext)
  plot(mybgxs, col=mycols1, box=FALSE, bty="n", cex.axis=0.8, legend=FALSE, axes=FALSE)
  plot(worldmap[5], border="grey20", lwd= 0.5, add=TRUE,col=NA,reset=FALSE)
  
  for (mapunit in SppSum[ns,"gen_sp_subtax"]){
    i<-which(SppSum[ns,"gen_sp_subtax"]==mapunit)
    k<-which(SppSum$gen_sp_subtax==mapunit)
    yesCountry<-str_split(SppSum[k,"country"], "/")                   #split countries of occurrence into separate elements
    mycountries <- subset(worldmap, worldmap$GU_A3 %in% yesCountry[[1]])     #get subset of countries that species occurs in
    plot(mycountries[5], add=TRUE, col=NA, border="black",lwd= 1)
    
    species_shapeA <- subset(WHOpolies, WHOpolies$NAME_SSN == gsub("_"," ",paste(mapunit)))
    #species_shapeB <- subset(WHOpolies, WHOpolies$NAME_SSN == gsub("_"," ",gsub(" 01","",paste(SppSum$WHOSpecies[k]))))
    species_shapeE <- subset(WHOpoliesNEW, WHOpoliesNEW$Scientific == gsub("_"," ",paste(mapunit)))
    #species_shapeF <- subset(WHOpoliesNEW, WHOpoliesNEW$Scientific == gsub("_"," ",gsub(" 01","",paste(SppSum$WHOSpecies[k]))))
    
    #if(!(length(species_shapeE$OBJECTID)>=1) & !(length(species_shapeF$OBJECTID)>=1)){ #if there's no updated polygons, plot old ones
    if(!(length(species_shapeE$OBJECTID)>=1)){ #if there's no updated polygons, plot old ones
      
      if(length(species_shapeA$OBJECTID)>=1){
        plot(species_shapeA[2], add=TRUE, border=mycols2[i],col=WHOcol,lwd= 1)
      }#else{
      #  if(length(species_shapeB$OBJECTID)>=1){
      #    plot(species_shapeB[2], add=TRUE, border=mycols2[i],col=WHOcol,lwd= 1)
      #  }
      #}
    } else { #otherwise use new versions
      
      if(length(species_shapeE$OBJECTID)>=1){
        plot(species_shapeE[2], add=TRUE, border=mycols2[i], col=WHOcol, lwd= 1)
      }#else{
      #if(length(species_shapeF$OBJECTID)>=1){
      #  plot(species_shapeF[2], add=TRUE, border=mycols2[i], col=WHOcol, lwd= 1)
      # }
      #}
    }}
  
  for (mapunit in SppSum[ns,"gen_sp_subtax"]){
    
    i<-which(SppSum[ns,"gen_sp_subtax"]==mapunit)
    if(file.exists(paste(MyDir, "/Maxent_SWDs/spp_occs/",mapunit,".csv",sep=""))){
      
      redrecs<-read.table(paste(MyDir, "/Maxent_SWDs/spp_occs/",mapunit,"_redRecs.csv",sep=""), sep=",", header=TRUE)
      points(redrecs$lat~redrecs$long, col=mycols2[i], pch=20, cex=1)
      points(redrecs$lat[redrecs$VetCom==2]~redrecs$long[redrecs$VetCom==2], col=mycols2[i], pch=4, cex=1)
      
      #add a species label to the centre of the record cloud for each mapping unit
      myx<-mean(redrecs$long)
      myy<-mean(redrecs$lat)
      graphics::text(x=myx, y=myy, labels = gsub("_"," ",mapunit), cex = 1, col = mycols2[i])
      
    }}
  
  #plot title in top margin of page:
  mtext(paste("Antivenom & human impact group:"), font= 2, side = 3, line = 1, outer = TRUE, cex= 1.4)
  mtext(paste("'",mapname,"'",sep=""), font= 4, side = 3, line = -0.3, outer = TRUE, cex= 1.4)
  mtext(paste("includes species:"), font= 1, side = 3, line = -1.5, outer = TRUE, cex= 1)
  
  if(length(ns)<=4){
    mtext(paste(SppSum[ns,"Accepted_Species"], collapse = ', '), 
          font= 3, side = 3, line = -2.4, outer = TRUE, cex= 0.9)
  } else {
    if(length(ns)<=8){
      mtext(paste(SppSum[ns[1:4],"Accepted_Species"], collapse = ', '), 
            font= 3, side = 3, line = -2.4, outer = TRUE, cex= 0.9)
      mtext(paste(SppSum[ns[5:length(ns)],"Accepted_Species"], collapse = ', '), 
            font= 3, side = 3, line = -3.5, outer = TRUE, cex= 0.9)
    } else {
      if(length(ns)<=12){
        mtext(paste(SppSum[ns[1:4],"Accepted_Species"], collapse = ', '), 
              font= 3, side = 3, line = -2.4, outer = TRUE, cex= 0.9)
        mtext(paste(SppSum[ns[5:8],"Accepted_Species"], collapse = ', '), 
              font= 3, side = 3, line = -3.3, outer = TRUE, cex= 0.9)
        mtext(paste(SppSum[ns[9:length(ns)],"Accepted_Species"], collapse = ', '), 
              font= 3, side = 3, line = -4.2, outer = TRUE, cex= 0.9)
      } else {
        mtext(paste(SppSum[ns[1:4],"Accepted_Species"], collapse = ', '), 
              font= 3, side = 3, line = -2.4, outer = TRUE, cex= 0.9)
        mtext(paste(SppSum[ns[5:8],"Accepted_Species"], collapse = ', '), 
              font= 3, side = 3, line = -3.3, outer = TRUE, cex= 0.9)
        mtext(paste(SppSum[ns[9:12],"Accepted_Species"], collapse = ', '), 
              font= 3, side = 3, line = -4.2, outer = TRUE, cex= 0.9)
        mtext(paste(SppSum[ns[13:length(ns)],"Accepted_Species"], collapse = ', '), 
              font= 3, side = 3, line = -5.1, outer = TRUE, cex= 0.9)
      }}}
  #fonts: 1=normal, 2=bold, 3=cursive, 4= bold & cursive
  
  dev.off()
  
}
}







############################
# make fourth lot of overview plots for species listing groups
############################

################
################
mapspp<-unique(SppSum$Listing_Name)
mapspp<-gsub(" ","_",gsub("; ","_",mapspp))
mapsppB<-unique(SppSum$Listing_Name[SppSum$Redo==1])
#mapsppB<-c("Walterinnesia aegyptia","Walterinnesia morgani")
#mapsppB<-unique(SppSum$Listing_Name)
mapsppB<-gsub(" ","_",gsub("; ","_",mapsppB))
################
################

#SppSum$maprecs<-paste("NA")

for (mapsp in mapspp){
  if(mapsp %in% mapsppB){
  
  ns<-which(SppSum$Listing_Name==gsub("_"," ",mapsp))
  m<-which(mapspp==mapsp)
  maprecs<-data.frame()
  mapname<-unique(SppSum$Listing_Name)[m]
  
  Sys.sleep(0.1)
  print(paste(Sys.time()," plotting data for WHO listing groupx ", gsub("_"," ",mapsp), ": ", m, " out of ", length(mapspp), sep=""))
  
  #make PDF plots:
  png(filename=paste(PlotDir,"/WHOListGroup_",mapsp,"_recs.png",sep=""), units="in", width=8.27, height=11.69, res=500) #png smaller file size, A3 for better resolution
  par(mfrow = c(2,1), mar = c(2.5,2.5,0,0), oma=c(1,1,4,3), mgp=c(1.2,0.5,0))
  
  plot(mybg, col=mycols1,box=FALSE, bty="n", cex.axis=0.8, legend=FALSE, axes=FALSE)
  plot(worldmap[5], border="grey20", lwd= 0.5, add=TRUE,col=NA,reset=FALSE)
  
  for (mapunit in SppSum[ns,"gen_sp_subtax"]){
    i<-which(SppSum[ns,"gen_sp_subtax"]==mapunit)
    k<-which(SppSum$gen_sp_subtax==mapunit)
    yesCountry<-str_split(SppSum[k,"country"], "/")                   #split countries of occurrence into separate elements
    mycountries <- subset(worldmap, worldmap$GU_A3 %in% yesCountry[[1]])     #get subset of countries that species occurs in
    plot(mycountries[5], add=TRUE, col=NA, border="black",lwd= 1)
  }
  
  #bind all records for species included in this map species together
  for (mapunit in SppSum[ns,"gen_sp_subtax"]){
    i<-which(SppSum[ns,"gen_sp_subtax"]==mapunit)
    
    if(file.exists(paste(MyDir, "/Maxent_SWDs/spp_occs/",mapunit,".csv",sep=""))){
      addrecs<-read.table(paste(MyDir, "/Maxent_SWDs/spp_occs/",mapunit,".csv",sep=""), sep=",", header=TRUE)
      maprecs<-rbind(maprecs,addrecs)
      
      redrecs<-read.table(paste(MyDir, "/Maxent_SWDs/spp_occs/",mapunit,"_redRecs.csv",sep=""), sep=",", header=TRUE)
      points(redrecs$lat~redrecs$long, col=mycols2[i], pch=20, cex=1)
      points(redrecs$lat[redrecs$VetCom==2]~redrecs$long[redrecs$VetCom==2], col=mycols2[i], pch=4, cex=1)
      
      #add a species label to the centre of the record cloud for each mapping unit
      myx<-mean(redrecs$long)
      myy<-mean(redrecs$lat)
      #graphics::text(x=myx, y=myy, labels = gsub("_"," ",mapunit), cex = 1, col = mycols2[i])
      
    }}
  
  if(length(maprecs)>0){
    l<-length(maprecs[,1])
    #SppSum$maprecs[ns]<-paste(as.character(l))
    
    #write out copy of reduced species records:
    maprecs<-subset(maprecs, !is.na(maprecs$myTaxon))
    #write.table(x=maprecs, file=paste(MyDir, "/Maxent_SWDs/spp_occs/",mapsp,"_mod.csv",sep=""),qmethod = "double",
    #            na = "NA", row.names = FALSE, col.names = TRUE, sep=",",quote=TRUE)
    
    #make extent object describing a box around occurrence records (+-5degrees):
    myext<-extent(min(na.omit(maprecs$long))-5,  
                  max(na.omit(maprecs$long))+5, 
                  min(na.omit(maprecs$lat))-5, 
                  max(na.omit(maprecs$lat))+5)
    
  } else {
    myext<-extent(mybg)
  }
  
  mybgxs<-crop(mybg, y=myext)
  plot(mybgxs, col=mycols1, box=FALSE, bty="n", cex.axis=0.8, legend=FALSE, axes=FALSE)
  plot(worldmap[5], border="grey20", lwd= 0.5, add=TRUE,col=NA,reset=FALSE)
  
  for (mapunit in SppSum[ns,"gen_sp_subtax"]){
    i<-which(SppSum[ns,"gen_sp_subtax"]==mapunit)
    k<-which(SppSum$gen_sp_subtax==mapunit)
    yesCountry<-str_split(SppSum[k,"country"], "/")                   #split countries of occurrence into separate elements
    mycountries <- subset(worldmap, worldmap$GU_A3 %in% yesCountry[[1]])     #get subset of countries that species occurs in
    plot(mycountries[5], add=TRUE, col=NA, border="black",lwd= 1)
    
    species_shapeA <- subset(WHOpolies, WHOpolies$NAME_SSN == gsub("_"," ",paste(mapunit)))
    #species_shapeB <- subset(WHOpolies, WHOpolies$NAME_SSN == gsub("_"," ",gsub(" 01","",paste(SppSum$WHOSpecies[k]))))
    species_shapeE <- subset(WHOpoliesNEW, WHOpoliesNEW$Scientific == gsub("_"," ",paste(mapunit)))
    #species_shapeF <- subset(WHOpoliesNEW, WHOpoliesNEW$Scientific == gsub("_"," ",gsub(" 01","",paste(SppSum$WHOSpecies[k]))))
    
    #if(!(length(species_shapeE$OBJECTID)>=1) & !(length(species_shapeF$OBJECTID)>=1)){ #if there's no updated polygons, plot old ones
    if(!(length(species_shapeE$OBJECTID)>=1)){ #if there's no updated polygons, plot old ones
      
      if(length(species_shapeA$OBJECTID)>=1){
        plot(species_shapeA[2], add=TRUE, border=mycols2[i],col=WHOcol,lwd= 1)
      }#else{
      #  if(length(species_shapeB$OBJECTID)>=1){
      #    plot(species_shapeB[2], add=TRUE, border=mycols2[i],col=WHOcol,lwd= 1)
      #  }
      #}
    } else { #otherwise use new versions
      
      if(length(species_shapeE$OBJECTID)>=1){
        plot(species_shapeE[2], add=TRUE, border=mycols2[i], col=WHOcol, lwd= 1)
      }#else{
      #if(length(species_shapeF$OBJECTID)>=1){
      #  plot(species_shapeF[2], add=TRUE, border=mycols2[i], col=WHOcol, lwd= 1)
      # }
      #}
    }}
  
  for (mapunit in SppSum[ns,"gen_sp_subtax"]){
    
    i<-which(SppSum[ns,"gen_sp_subtax"]==mapunit)
    if(file.exists(paste(MyDir, "/Maxent_SWDs/spp_occs/",mapunit,".csv",sep=""))){
      
      redrecs<-read.table(paste(MyDir, "/Maxent_SWDs/spp_occs/",mapunit,"_redRecs.csv",sep=""), sep=",", header=TRUE)
      points(redrecs$lat~redrecs$long, col=mycols2[i], pch=20, cex=1)
      points(redrecs$lat[redrecs$VetCom==2]~redrecs$long[redrecs$VetCom==2], col=mycols2[i], pch=4, cex=1)
      
      #add a species label to the centre of the record cloud for each mapping unit
      myx<-mean(redrecs$long)
      myy<-mean(redrecs$lat)
      graphics::text(x=myx, y=myy, labels = gsub("_"," ",mapunit), cex = 1, col = mycols2[i])
      
    }}
  
  #plot title in top margin of page:
  mtext(paste("WHO Listing group:"), font= 2, side = 3, line = 1, outer = TRUE, cex= 1.4)
  mtext(paste("'",mapname,"'",sep=""), font= 4, side = 3, line = -0.3, outer = TRUE, cex= 1.4)
  mtext(paste("includes species:"), font= 1, side = 3, line = -1.5, outer = TRUE, cex= 1)
  
  if(length(ns)<=4){
    mtext(paste(SppSum[ns,"Accepted_Species"], collapse = ', '), 
          font= 3, side = 3, line = -2.4, outer = TRUE, cex= 0.9)
  } else {
    if(length(ns)<=8){
      mtext(paste(SppSum[ns[1:4],"Accepted_Species"], collapse = ', '), 
            font= 3, side = 3, line = -2.4, outer = TRUE, cex= 0.9)
      mtext(paste(SppSum[ns[5:length(ns)],"Accepted_Species"], collapse = ', '), 
            font= 3, side = 3, line = -3.5, outer = TRUE, cex= 0.9)
    } else {
      if(length(ns)<=12){
        mtext(paste(SppSum[ns[1:4],"Accepted_Species"], collapse = ', '), 
              font= 3, side = 3, line = -2.4, outer = TRUE, cex= 0.9)
        mtext(paste(SppSum[ns[5:8],"Accepted_Species"], collapse = ', '), 
              font= 3, side = 3, line = -3.3, outer = TRUE, cex= 0.9)
        mtext(paste(SppSum[ns[9:length(ns)],"Accepted_Species"], collapse = ', '), 
              font= 3, side = 3, line = -4.2, outer = TRUE, cex= 0.9)
      } else {
        mtext(paste(SppSum[ns[1:4],"Accepted_Species"], collapse = ', '), 
              font= 3, side = 3, line = -2.4, outer = TRUE, cex= 0.9)
        mtext(paste(SppSum[ns[5:8],"Accepted_Species"], collapse = ', '), 
              font= 3, side = 3, line = -3.3, outer = TRUE, cex= 0.9)
        mtext(paste(SppSum[ns[9:12],"Accepted_Species"], collapse = ', '), 
              font= 3, side = 3, line = -4.2, outer = TRUE, cex= 0.9)
        mtext(paste(SppSum[ns[13:length(ns)],"Accepted_Species"], collapse = ', '), 
              font= 3, side = 3, line = -5.1, outer = TRUE, cex= 0.9)
      }}}
  #fonts: 1=normal, 2=bold, 3=cursive, 4= bold & cursive
  
  dev.off()
  
}
}



##################################################
#write out copy of new species summary file:
write.table(x=SppSum, file=paste(MyDir, "/Maxent_LUTables/",projectName,"_Fin1b.csv",sep=""),
            row.names = FALSE, col.names = TRUE, sep=",")

############################################
############################################
# next: 3-VAP_MaxentPrep ###################
############################################